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Abstract. We show that the formation of large-scale structures through gravitational instability in the expanding 
universe can be fully described through a path-integral formalism. We derive the action S[f] which gives the 
statistical weight associated with any phase-space distribution function /(x, p, t). Fhis action S describes both 
the average over the Gaussian initial conditions and the Vlasov-Poisson dynamics. Next, applying a standard 
method borrowed from field theory we generalize our problem to an JV— field system and we look for an expansion 
over powers of 1/N. We describe three such methods and we derive the corresponding equations of motion at the 
lowest non-trivial order for the case of gravitational clustering. This yields a set of non-linear equations for the 
mean / and the two-point correlation G of the phase-space distribution /, as well as for the response function R. 
These systematic schemes match the usual perturbative expansion on quasi-linear scales but should also be able 
to treat the non-linear regime. Our approach can also be extended to non-Gaussian initial conditions and may 
£ — ■ serve as a basis for other tools borrowed from field theory. 
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q. 1. Introduction 

i • 

Q ■ In standard cosmological scenarios large-scale structures in the universe have formed through the growth of small initial 
density fluctuations by gravitational instability, see Peebles (1980) Once they enter the non-linear regime, high-density 
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^2 fluctuations build massive objects which will give birth to galaxies or clusters. Thus, it is important to understand 
. . \ the dynamics of these density perturbations in order to describe the large-scale structures we observe in the present 
universe. Moreover, observations (e.g. weak-lensing or galaxy surveys) can provide a measure of the power-spectrum 
of the density fluctuations down to non-linear scales. Hence it is useful to devise theoretical tools which may treat the 
±2 non-linear regime to make the connection between the data and the primordial power-spectrum. 

However, the non-linear Vlasov-Poisson equations that describe the dynamics of such density fluctuations within 
an expanding universe are not easy to solve and most rigorous approaches are based on perturbative expansions. 
More precisely, the distribution of matter is usually described as a pressure-less fluid through the standard equations 
of hydrodynamics (continuity and Euler equations) coupled to the Poisson equation for the gravitational potential. 
Then, one often looks for a solution of these equations of motion in the form of a perturbative expansion over powers 
of the initial fluctuations, or rather over powers of the linear growing mode (e.g., |Fry (1984)| |Goroff et al. (1986)1 



Bernardeau et al. (1994)1. An alternative method which is not based on such perturbative expansions nor on the 



hydrodynamical approximation, developed in Valageas (2002a'){ can go somewhat beyond these previous approaches 



for the statistics of the density within spherical cells. However, all these methods break down at shell-crossing and 
cannot handle the non-linear regime. Hence they are restricted to large scales, except for the approach described in 



Valageas (2002b) which can still handle very rare voids in the non- linear regime (where shell-crossing is not important). 



An alternative approach to this problem is to investigate simplified equations of motion which would capture most 



of the physics and would be easier to analyze. Thus, the well-known Zel'dovich approximation I Zel'dovich (1970) I 
extrapolates the trajectories of the particles from the linear regime. Next, the density field is reconstructed from the 
particles displacements. The advantage of looking at the trajectories is that one can handle large density perturba- 
tions and even go beyond shell-crossing: particles may cross each other and their trajectories remain well-defined while 
the local density momentarily diverges (which prevents the application of naive perturbative methods to the density 
field itself). Numerical studies have shown that the Zel'dovich approximation works well until shell-crossing and it 
can reproduce the large-scale skeleton of the matter distribution (e.g., Coles et al. (1993) I. However, it breaks down 
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after shell-crossing since it does not include any feedback: particles keep moving ahead after pancake formation and 
are not pulled back by the gravitational potential well. This has led to the adhesion model (Gurbatov et al. (1989) 



Vergassola et al. (1994) I where one adds an artificial viscosity term to the Euler equation so as to mimic the gravita- 
tional "sticking" of particles within potential wells. Then, the thickness of filaments or halos depends on this viscosity 
parameter. This model can be generalized by including a dependence on density and time within the viscosity term 



(Buchert et al. (1999) I which may improve the predictions obtained for the thickness of non-linear objects. However, 



it remains to be seen whether one can describe the inner structure of non-linear halos by following this line of attack. 



A different approach was proposed by Widrow & Kaiser (1993) who suggested to replace the Vlasov-Poisson sys- 



tem by a Schrodinger-Poisson system. This model is much closer to the original Vlasov-Poisson system than the 
hydrodynamical approaches recalled above and it contains more phase-space information. Indeed, from the wave- 
function ip(x,t) one may construct a full distribution function /(x, p, t) (which is not restricted to a Dirac or a 
Maxwellian in velocity space) which follows an equation of motion that reduces to the Vlasov equation in some limit 



(see Widrow & Kaiser (1993) I. The advantage of using i/;(x, t) is that it reduces the phase space from six coordinates 
(x, p) down to three (x). As described in Widrow & Kaiser (1993) and Davies fc Widrow (1997)[ this approach can 



treat shell-crossing and virialization. A comparison with the Zel'dovich approximation and simple examples show that 



this method is indeed a promising tool for cosmological purposes (Coles & Spencer (2003) I. 

In this article, we develop a new approach to this problem by working directly with the original Vlasov-Poisson 
equations. For cosmological purposes, we are not really interested in the evolution of the system for a few specific 
initial conditions but we rather wish to derive the statistical properties of the system after averaging over these initial 
conditions (usually taken to be Gaussian from inflationary scenarios). This suggests the use of a functional approach 
as in usual statistical physics, which directly focuses on the correlation functions of the system (or other statistics). 
However, in order to obtain such quantities one must resort to approximation schemes since usually non-Gaussian 
path-integrals cannot be computed exactly. Not to mix up the sources of error introduced by such approximations 
(associated with the procedure used to derive averaged quantities) with those entailed by the use of approximate 
equations of motion it is desirable to first apply this approach to the exact Vlasov-Poisson system. Moreover, the 
ultimate goal is to be able to describe the original system itself. However, our study is also complimentary to previous 
works which investigated hydrodynamical or "quantum" systems as a substitute to the Vlasov-Poisson equations since 
it could be applied to such approaches. 

Thus, we show that the statistical properties of the system, assuming Gaussian initial conditions, are fully described 
by a path-integral formalism and we derive the action S[f] which gives the statistical weight e~ s ^ associated with 
any phase-space distribution function /(x, p, t). This action S contains both the average over the Gaussian initial 
conditions and the Vlasov-Poisson dynamics. Moreover, we note that it yields the response functions of the system in 
addition to the correlation functions. Then, this result allows us to use standard tools of field theory. In particular, in 
this article we describe three such methods, based on "large- V expansions, and we derive for the case of gravitational 
clustering the equations of motion they imply for the mean phase-space distribution / and its two-point correlation 
G, at next-to- leading order (the first non-trivial order). 

This paper is organized as follows. In Sect we recall the Vlasov-Poisson equations of motion and we introduce our 
notations. We also pay some attention to the infrared divergences which appear for power-law linear power-spectra 
Ph{k) oc k n with n < — 1. Next, in Sect[2]we show how our system can be described through a path-integral formalism 
and we derive the corresponding action S[f]. We also discuss the response functions of the system and we briefly 
review several approximation schemes. Then, in Sect0]we derive up to next-to-leading order the equations given by a 
direct steepest-descent method. We discuss in Sect03an alternative 1PI effective action scheme and we finally obtain in 
Sect the Schwinger-Dyson equations given by a third approach based on the 2PI effective action (we also present in 
App[X]a simple proof of a key result needed in this approach). We also show that the infrared behaviour encountered 
in the case of large-scale structure formation selects some approaches over others. Finally, we investigate a standard 
perturbative expansion in Sect d and we discuss our results in Sect|Hl 

2. Equations of motion 

2.1. Vlasov-Poisson system 

First, we recall here the equations of motion which describe the gravitational dynamics of a collisionlcss fluid in an 
expanding universe. We consider homogeneous initial conditions, that is the system is statistically invariant through 
translations. In the continuous limit where the mass m of the particles goes to zero the system is fully described by 
the distribution function (phase-space density) /(x, p,i) where we note x the comoving coordinate of a particle and 
p its momentum (which we shall also refer to as the velocity) defined by: 



a 2 x, 



and we normalize / by: p(x, t) = p J /(x, p, t) dp. (1) 
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Here we introduced the scale factor a(t), the comoving density p(x, t) and the mean comoving density of the universe 
p (which is independent of time). Note that we normalized the distribution function /(x, p, t) by the constant mass 
density p, so that J /dp is dimensionless. Then, the gravitational dynamics is given by the collisionless Boltzmann 
equation, also called Vlasov equation, coupled with the Poisson equation ( Peebles (1980) ) : 



df p df d^ df AnG AirG p(x, t) - p 

— + ^.- — .— =0, with A$ = (p- p ) = p 5 and <5(x, t) = , (2) 

at a z ox ox op a a p 

where $ is the gravitational potential and S is the density contrast. The fact that the Poisson equation involves 
the density perturbation (p — p) rather than the total density p comes from the change from physical coordinates 
to comoving coordinates (this also removes the cosmological constant term Q\ from these equations of motion). In 
order to derive eq.© we neglected the discrete character of the distribution of matter and we assumed small-scale 



smoothness (i.e. there is a small-scale cutoff to the power-spectrum of the density fluctuations), see Peebles (1980) 
As can be seen from eq.lPJ, the gravitational force F(x, t) can be written as: 

F(x,t) = = g /dx' (p(x')-p) ^— * = % [ dx'-p8{x>) ^— * = g /dx' p(x') J^ZJL. (3) 



dx 

The fourth term in eq.J3J) merely uses the definition of the density contrast S(x, t) while the fifth term uses the fact 



that the integration over the angles of (x' — x) of the factor multiplied by the constant ~p vanishes (see Peebles (1980) I. 
Note that this actually involves an implicit temporary regularization of the gravitational interaction at large scales (see 
eq.JfJJl below). This is sometimes referred to as the "Jeans swindle". For practical purposes, it is convenient at some 
points to work in global Fourier space. Thus, we introduce the Fourier conjugates k and w of the spatial coordinate x 
and of the momentum p, and we define the Fourier transforms of the distribution function f(x, p, t) and of the density 
contrast 5(x,t) by: 

/(x,p,t) = y dkdw e i(k x+w p) /(k,w,<) and 5{x, t) = J dk e ik x S(k, t). (4) 

To avoid introducing too many variable names, we use the same notations both in real space (x, p) and in Fourier 
space (k, w). The arguments should remove any ambiguity. We also define the Fourier transforms of other fields as in 
eq.J3J. In Fourier space, the gravitational force can be obtained from eq.© or eq.Q and we write: 

F(k, t) = i^£p(k, t) = l ^(2^) 3 ^/(k, 0, t). (5) 

Note that we could have chosen the density perturbation {p — ~p) rather than the total density p in eq.(|5j), as can be 
seen from eq.Q. This merely translates the fact that the gravitational potential is not fully defined by the Poisson 
equation (one should add boundary conditions). This is directly related to the "Jeans swindle" discussed below eq.(0l 
which also leads to the "property" : 

k k k 

— ~ Sn (k) = 0, because —5- stands for lim —5 -5-, (6) 

K l k l k a ^0 K l + hp 

where 80 stands for Dirac's distribution and we wrote explicitly the large-scale regularization of the gravitational 
interaction (beyond the scale l/k c — > 00). In practice we do not need to write explicitly the cutoff k c and all our 
equations are written with k c = (i.e. no large-scale cutoff) because we usually do not encounter ambiguous quantities 
as in JBJ. The only exceptions are ea. H19(l and ea.l|48|) below where we must remember the "property" Then, using 
eq.JSJ the Vlasov equation J2J) reads in Fourier space: 

df_kd£ = 4^p (27r)3 f dk , k^w 
dt a 2 dw a J k 

This non-linear and non-local equation describes the dynamics of the collisionless fluid. Thus, in order to investigate 
the formation of large-scale structures in the universe (neglecting the influence of baryons onto the dark matter) we 
could solve eq.lQ) for some specific initial conditions and next perform a suitable average over these initial conditions 
(characterized by the power-spectrum of the early linear perturbations). In the following, introducing the matrix O 
and the vertex K , we shall write the quadratic equation (JJJ) in the more concise form : 

0{xi,x 2 )f(x 2 ) = K(xi;x 2l x 3 )f(x 2 )f{x 3 ) + 7ft (xi), with f(xi) = 6(ti - U)f(xi) and rfr{xi) = 6 D (ti - ii^jfei), (8) 

where the argument x = (k, w, t) = (x, t) (or x = (x, p, t) in real space) represents the space, velocity and time 
variables, the set of all variables excluding time is denoted by underlining the argument (e.g., x); and we used the 
integration convention for repeated indices. The new term iji(xi) in the r.h.s. in eq.JSJ) describes the initial conditions 
at time tj so that f(x) = for t < t; and f(x, t;) = n.(x). We shall eventually take the limit U — > 0. Throughout this 
article we note Sd and 6 the Dirac and Heaviside functions. 



(2tt) 3 /dk'— r /(k',0)/(k-k',w). (7) 
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2.2. Initial conditions 



In this paper we consider standard cosmological scenarios where large-scale structures have formed through the growth 
of small density fluctuations by gravitational instability. Therefore, at early times the system shows small perturbations 
from the homogeneous expanding universe and one can use linear theory. In practice, this regime is usually studied 
within the framework of the hydrodynamical approximation where all particles at a given point x have the same 
momentum p(x) and the linear growing mode is given in partial Fourier space by (Peebles (1980)): 



<5 L (k,i) = D+(t)5 L oQi) and p L (k,i) = ia 2 D+ — 5 L o(k), with D 



2 a -D + = 



AttQp 



(9) 



Here 6l and are the linear density contrast and momentum while i5lo( x ) is the linear density contrast today (at 
z = 0) and D + (t) is the linear growth factor. Next, one can look for a solution to the Euler-Poisson system through 
a perturbative expansion over powers of the linear growing mode Sl (e.g., Bernardeau et al. (2002) for a review). As 
seen in Valageas (2001) the same procedure can be applied to the Vlasov-Poisson system and it yields identical results 



at all orders (which strongly suggests that shell-crossing is beyond the reach of such perturbative schemes). The linear 
growing mode t)l of the Vlasov equation which corresponds to © is given by (see Valageas (2001) I: 



77i(x,p,t) = 5 D (p) + 6 L (x,t)5 D (p) 



-PL{x,t). — {p) and r) L {k,w,t) = -j—^ 



(2tt) 3 



<M k ) 



k.w 



(10) 



Here we included the constant term <5d(p) which corresponds to the homogeneous expanding universe. Thus, at 
early times (or large scales) we have / — ► t)l- The Dirac factors <5d(p) in ea. (|10fl translate the fact that we have no 
velocity dispersion at linear order (we have a pressure-less hydrodynamical fluid) and higher orders of the perturbative 
expansion over powers of the linear mode t]l merely generate higher-order derivatives of Sn(p) (see Valageas (2001) I. 
Then, we can choose for the initial conditions Tfi introduced in eq.©: 



77.(x, p) = ?/ L (x,p,ii), and we note = <r?i> and Ai(xi, x 2 ) = {rj- 1 (x 1 )rj i (x 2 )) c , 



(11) 



where (..} denotes the average over the initial conditions. In this article we shall consider Gaussian initial conditions, 
so that both the linear density field 6l and the linear growing mode r\L are Gaussian: 



(<W k i)<Wk 2 )} = ^Lo(&i)<M k i +k 2 ) and we note A L (xi,x 2 ) = {vl(xi)ti l (x 2 )} c . 



(12) 



Then, the field rji is also Gaussian and its connected two-point correlation A ; can be expressed in terms of the linear 
power-spectrum Pto(fc) through ea. (|10f) . Note that for a non-zero initial time U > the choice (|11H for the initial 
conditions r\. actually yields both growing and decaying linear modes. However, in the limit t\ — ► we indeed recover 
the pure linear growing mode tjl defined in eg. 1101) . 



2.3. Long-wavelength divergences 

From eq.@, one can see that for a power-law linear power-spectrum Plo(^) °c k n (with —3 < n < 1 in the cosmological 
context) the rms linear momentum at scale R formally scales as tx p (i?) oc R~( n+1 ^ 2 . This actually leads to infrared 
divergences for the linear velocity field p^ as defined in eq.@ for n < — 1. More precisely, if we try to set up some 
generic initial conditions from eq.JUJl we can see that long wavelengths yield a divergent contribution to p^(x). Thus, 
we can no longer use eq.© to define the linear velocity field. Within the standard perturbative approach based on 
the Euler-Poisson system, these infrared divergences actually cancel out (at least at leading order) when one only 
considers the equal-time correlation of the density field, as seen in Jain & Bertschinger (1996) (also Vishniac (1983) I. 



As pointed out by these authors, this can be understood from the fact that contributions to the velocity field from 
long-wave modes correspond to an almost uniform translation of the fluid and therefore should not affect the growth 
of density structures on small scales. However, these infrared contributions can no longer be ignored when one studies 
different-times correlations or the velocity field itself, which is clearly the case when we work with the Vlasov equation. 
As suggested in Jain &: Bertschinger (1996)1 the only quantity which really matters is the velocity difference between 



neighbouring points and not the overall mean velocity. This is not surprising since the equations of motion written 
in real physical space obey the usual Galilean symmetry (note that this is no longer the case for eq.(0) or eq.JTJ) 
which are written in comoving coordinates). Therefore, we must recover the same physical behaviour if we subtract 
the initial velocity field by some constant. In particular, we could normalize the linear velocity so that p^(x = 0) = 
for instance. Thus, we now define the linear velocity field Pl(x) by: 

Pl(x, t) = fdk (e ik x - 1) ia 2 D+ ± S L0 (k). (13) 
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One can sec that this still yields a linear solution to the Vlasov equation provided we subtract a constant term over 
space to the gravitational force so that F(x = 0) = 0. Of course, this simple change of normalization does not modify 
the physics of gravitational clustering and it is consistent with Poisson's equation. Finally, this yields for the Vlasov 
equation: 

ft-i^ = 4j ?^I dk ' W<- - ^ • < 14 > 

By comparison with eq.Q) we see that the "counter-term" /(k, w) in the bracket ensures that the contribution from 
large scales (i.e. k' — > 0) converges (when we work in real space x). Therefore, ea. l|14|) allows us to investigate linear 
power-spectra with n < — 1. Let us stress that ea. <|14|) is exact and is strictly equivalent to eq.(|7J) with respect to 
gravitational clustering, that is the density fields obtained from both evolution equations are identical. However, the 
velocities and the positions of the particles differ between both frameworks by some values which are constant over 
space (but change with time). These offsets actually diverge for n < — 1 because of long wavelength contributions but 
they have no influence on the physics at work. Working with ea. (|14f> ensures that when we investigate the physics at 
some scale R we only encounter quantities which are governed by this scale (or possibly the scale which is turning 
non- linear), and not by some much larger scale (where n crosses the threshold —1) which eventually disappears as 
various contributions cancel out. This is a priori of great interest for numerical purposes. Unfortunately, within this 
framework we lost homogeneity (i.e. invariance through translations) since we singled out the origin x = 0. This 
makes the analysis of the problem and the interpretation itself of the correlations of the density field more intricate. 
Therefore, in the following we shall only consider the case of homogeneous initial conditions (i.e. with a large-scale 
cutoff) which should be sufficient for practical purposes. However, the influence of these infrared divergences will play 
a key role in Sect 16.41 as they will discriminate between various approximation schemes. 



3. Functional formulation 

3.1. Derivation of the actions S[f] and S[f, A], mean response functions 

In practice, we are not really interested in the exact solution / of the equation of motion for a given initial condition rji . 
Indeed, since the initial condition rji is a random field we are mainly interested in the statistical properties of the phase- 
space distribution function /. They can be described through a Heisenberg operator formalism (Martin et al. (1973) I 
or a functional approach (e.g., Phythian (1977)} Jensen (1981) I. We shall apply to the Vlasov-Poisson equation I 
the path-integral formalism which is more convenient. Hereafter we work in real space (x, p, t) (except in eas. ()17t - (fTH I 
below and in Sect[7J|. Thus, let us consider any functional F(f) of the classical field / (e.g., an n— point correlation). 
Its average over the initial conditions can be computed as: 

(F(f)) = [ [d m ] F(f) e-^-^A^-ta-n,) = f [d/] | DetM | e -h(Of-K a f-rj i ).ArK(of-K a f>-v i ) ) 



(15) 

where we have taken the Gaussian average over r/i and in the last term we made the change of variable rji — > f, 
using eq.JSJ. We also introduced the symmetric vertex (over its last two arguments) K s (x; X\, x%) = [K(x;xi,x 2 ) + 
K(x;x2,x±)]/2. We do not write explicitly the normalization factors which appear in front of path-integrals like (|15|) . 
If needed, they can be recovered from the constraint (1) = 1. The Jacobian |DetM|, with M — Srji/Sf, is given by: 



DetM = Det(C - 2KJ) = Det(C)Det(l - 2KJ) = Det(O) exp[Trln(l - 2KJ)], 



(16) 



where 1 — 2K s f is the functional derivative with respect to / of the integral form of the Vlasov equation Q)-© (see 
also |Zinn- Justin (1989)] )- The latter can be integrated by the method of characteristics which yields: 



/(k,w,t) = 77. [k,w 

+ 47rgp(27r) 3 



K ' a' 2 

4 oT 



dk' 



k' 2 ' 



dt^_ 
^2 



/(k',0,t')/ k-k',w + k 



d^ 
~d j2 '' 



so that the vertex K reads in Fourier space (with all times larger than t\): 



K(x;x 1 ,x 2 ) 



(2tt)^ - h) S D (h - t 2 ) fo(w x ) S D w 2 - w - k 



7J 2 



ki-w 2 



fc(ki +k 2 - k). 



(17) 



(18) 



Of course, we have O.K. = K, and we defined the symmetric vertex K s (x; x\, x 2 ) = 
Next, we expand the logarithm in ea. p6|l which yields: 

00 „ 
Trln(l - 2K s f) = - ^ - Tr (2K s f) q = -Tr 2K s f = since / dx K s (x; x,y) = 

q=l q J 



[K(x;xi,x 2 ) + K(x;x 2 ,Xi)]/2. 



(19) 
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Here we used the fact that Tr(2K s f) q = for q > 2 (see also Valageas (2001) Zinn- Justin (1989) I because of causality 
(i.e. the Heaviside factor which appears in the kernel K s ) and the last equality uses the property ||BJ) to remove any 
ambiguity. Therefore, the Jacobian in ea. l|15|l is equal to unity, up to the irrelevant constant Det(O). Thus, we obtain: 

(F(f)) = J [d/] F(f) e- s ^ with S[f] = \{Of- K s f - %) .A." 1 . {Of - K s f r},) . (20) 

Thus, ea. l|20|l shows that our system is fully described by the action S[f] therein defined. Note that S[f] depends 
on the configuration of the field / over phase-space (x, p) and time t. In particular, it is clear that our problem is 
non-stationary: we study an out-of-equilibrium dynamics driven by gravitational instability. Note that the procedure 
described above could easily be extended to non-Gaussian initial conditions, provided the statistical weight of the field 
r]i is still given by an action S^J^i], which would no longer be quadratic. From eas. (|10fl - H12|) one can easily show that 
the matrix A; is positive. Therefore, its inverse Ar 1 is also positive, hence the path-integral (|20|l is well-defined since 
S[f] > (the action S admits a lower bound). 

As we shall see below, it is convenient to write the path-integral 120|) in the form: 

(F(f)} = J[df][dX]F(f)e- s ^^ with S[f, A] = X.(Of - K s f 2 — rj^) - ^A.Aj.A, (21) 

where we introduced the auxiliary imaginary field A (i.e. A(a;) = iX'(x) with A' real). Indeed, we can check that the 
Gaussian integration over A of ea. (|21|l yields back cq. I|2U|I . One interest of the formulation l|21l) is that the auxiliary field 



A generates the response functions of the system ( |Phythian (1977) 1. Indeed, let us add a small non-random perturbation 
C(xi) to the right-hand side of the Vlasov equation Q)-||SJ). This amounts to the change 771 — > Tfi + C which can be 
described by the simple perturbation rji + C- Therefore, in ea. (|21|l the action is changed to S[f, X] — > S[f, A] — X.(. 
Then, the functional derivatives with respect to £ at the point £ = write: 

' W (/) -) C =o = / [d/][dA] X(x) F(f) e~ s ^ x \ <_^l_) c=0 = / [ d /][dA] X(x 1 )X(x 2 ) F(f) e' 3 ^,.. (22) 



This clearly shows that insertions of the auxiliary field A yield the mean response functions of the system. In particular, 
the mean first-order response function R(xi,x 2 ) is: 



R(xi,x 2 ) ~ 



= (f(x 1 )X(x 2 )) cx 6(h - i a ) while (X(x)) = 0, (A(xi)A(a; 2 )) = 0. (23) 



The last two equalities in (|23|l are obtained with F = 1. The response function R(x\,X2) is proportional to the 
Heaviside factor 9(ti — 1 2 ) because of causality: the field /(21) at time t\ only depends on the "noise" r}i{x 2 ) at earlier 
times t 2 < t\. Besides, the response R obeys the constraints: 

for ti — >■ tf : R(xi,x 2 ) — > Soi^ — x 2 ), and for all t\ > t 2 : J dxidpi R(xi,x 2 ) = 1, (24) 

where the last property translates the fact that J dxdp / is a constant of motion of the Vlasov equation. Since our 
system is an out-of-equilibrium dynamics there is no fluctuation-dissipation theorem: there is no linear relationship 
between the response R and the correlation G{x\,x 2 ) = {f(xi)f(x 2 )) c . 



3.2. Approximation schemes 

The path-integrals |(2Cj|| and (|2~TJ> provide a convenient approach to investigate the statistical properties of the system. 
However, usually one does not know how to compute such non-Gaussian path-integrals so that one must resort to 
various approximation schemes. A first approach is to apply a perturbative expansion by expanding the exponential 
of the non-Gaussian part of the action (this would give a series over powers of K s ). Such a perturbative scheme was 
described in Valageas (2001)[ based on an integral form of the Vlasov equation which directly relates 5f to the linear 
growing mode tjl (and not to the initial conditions at a finite time ti), where Sf = f — <$d(p) is the deviation of the 
phase-space density / from the homogeneous Hubble flow. This yields another action S[5f] whose non-Gaussian terms 
may be expanded. However, as shown in Valageas (2001)| this procedure actually recovers the standard perturbative 
results derived from the hydrodynamical approach. Therefore, it is not sufficient for our purpose which is to devise a 
theoretical tool to investigate the non-linear regime. On the other hand, a non-perturbative scheme is provided by the 
Feynman variational method (e.g., Feynman (1972) \ based on the inequality: 



[df] e- s[f] > e - {s - So)o / [df] e- So[f] with (S - 5*o)o 



J[df] (S[f]-S [f})e- S Q^ 
J[df] e~ s °if] 



(25) 
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This applies to any real actions S and So. Then, by maximizing the right-hand side of the inequality l|25[l over 
a class of Gaussian actions Sq where all integrals can be computed one may obtain a useful approximation. One 
may try to apply this method to the real action S[f] derived in ea. l(2U|l (note that S[f, A] in ea. (|2*T|l is complex). 
Unfortunately, this method cannot be easily applied to our system because the matrix A^ 1 is actually divergent since 
Aj is not invertible, as can be seen from eas. Hl(J|) - (|llfl . Thus the matrix A. is actually regularized in the action S[f] 
through an infinitesimal perturbation and we must make sure that such divergences disappear in the final results. 
This behaviour merely translates the fact that the random field r/i is constrained to be of the form This also 

means that the path-integral (|20|l only runs over phase-space distributions /(x, p, t) which are solutions of the Vlasov 
equation Q with initial conditions of the form Thus, in order to keep the quantity (S — So)o finite in 1)25(1 we 
must only consider trial actions Sq which also select solutions to the Vlasov equation with similar initial conditions. 
Therefore, we have not gained much by using this method: using such an action So is usually just as difficult as the 
original problem with the action S. Nevertheless, one might think to use a trial action Sq which only "authorizes" 
spherical linear density fields. Indeed, the gravitational dynamics of spherical fields is more easily solved (or at least 
easier to investigate) so that one could hope to derive in this way some useful information. However, although the 
average {S — £0)0 is well-behaved (since we integrate over exact solutions of the Vlasov equation) the path- integral 
in the right-hand side of the inequality I123H now involves a divergent Determinant. This comes from the fact that 
by integrating over spherical fields only we no longer have the same number of degrees of freedom in the systems 
defined by S or Sq- Therefore, we cannot derive any rigorous information from the restriction of the original problem 
to spherical fields only (see also Sect. 4 and App.A in Valageas (200 2b)| ). 

A third approach to path-integrals is to apply a steepest-descent expansion. Therefore, we may add a multiplicative 
parameter N to the action S and define the generating functional Zn [j, h] : 



Z N \j,h]= [df][dX 



^N[j.f+h.\-S(f,\)] 



whence Z N [j] = / [d/] e 



:; W[j./-I(0/-K s / 2 -7T i ).Ar 1 .(0/-/f a / 2 -7T 1 )] 



(26) 



As is well-known, functional derivatives of Zjq with respect to the source fields j and h yield all many-body correlations. 
The actual case corresponds to N = 1 but we may compute Zpj through an expansion over powers oll/N and next 
put N = 1 into the results. This procedure is similar to the "semiclassical" (or loopwise) expansion over powers of h 
encountered in quantum field theory. Note that this scheme is very different from the standard perturbative expansion 
described above since we make no difference between the quadratic and higher-order terms of the action, as the factor 
N multiplies the whole action S. 

Finally, a fourth method often used in physics to investigate such systems is to generalize our problem to an 
N— field system and look again for an expansion over powers of 1/N. Of course, there are many ways to generalize the 
action S[f] to N fields. In order to obtain good scaling properties with N, one can check that it is convenient to first 
write the partition function Z defined by the action S[f] derived in (|2(J[I as: 



Z= / [d/][dA][da] 



,-i(0/-77 i ).Ar 1 .(0/-7T i ) + A.(K s / 2 -<r) + (0/-r, i ).A i - 



r-icr.A.r 



(27) 



where we introduced the auxiliary fields A and er, with A imaginary. Next, we generalize our system from one field / 
to N similar fields /„ with 1 < a < N. Thus, we write the partition function Z' N as: 

r N 

Z' N = f[[df a }[d\}[da]e-iZj° f «-^- A "< ^ (28) 

J a=l 
, N 

= / I]j d ^ndA] e"^J /»^)- A r 1 -(o/o-^)+A.E a x a / a 2 +i[Ar 1 -E a (°/»-^)- A ]- A '-[ A r 1 -E 6 ( ^-^)- A ], (29) 

^ a=l 

where we performed the Gaussian integration over the field a in the last term. Then, in order to ensure that all terms 
contribute to the same order over N to the action, with the fields f a and A of order unity, we add some factors N in 
the exponent and we define the action Sn by: 



S N [fa, \] = l - i/O-Af 1 W« ~ ^ 



.AT 1 . 



b - Vi) 



-A. £(0/ - K a f a - 7/0 - ^A.Ai.A 



(30) 



Of course, for N = 1 we recover the action which corresponds to the case of only one field /, which is the case of 
practical interest to us. One can easily show that the real part of the action Sat [/a, A] is again positive. Note that the 
first two terms of Sn cancel out for N = 1, which shows that there are many ways to generalize to N fields, some 
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of which being more appropriate. In fact, these two terms are necessary for the path-integral Zn to be well-defined. 
Otherwise, the action Sn would only depend on the two fields A and a = ^ a {Of a — K s f% — rj-/) which would lead to 
divergences and inconsistencies (the system would not have truly been generalized to an (N + 1)— field problem). Note 
that the action Sn scales as Sn ~ N for large N (with all fields of order unity, compare also with ea. (|26f l). 

4. Direct steepest-descent method 

4.1. General formalism 

In the following, we shall describe how to apply to gravitational clustering three large- N methods which have been 
introduced in field theory. We shall consider the "semiclassical" ea. (|26[) which provides simpler expressions than the 
iV-field generalization i|30|) . The latter actually provides the same results (for N = 1) if we pay attention to the long- 
wavelength divergences discussed in Sect 12.31 see Sect 16.41 Although the action S[f, A] involves two fields we shall first 
recall the application of large- N expansions to the case of a single scalar field <f>(x) governed by a cubic action S[4>]: 

Z[j] = e NW ® = J [ d <p] e N ^~ s ^ with S[<P] = StixMx) + \s 2 {x, y)<f>(x)<P(y) + ^S 3 (x, y, z)^x)^y)<jy{z) (31) 

This corresponds to our case since the action S[f, A] is indeed cubic. This will yield simpler expressions in the inter- 
mediate steps and we shall directly translate the final results to our case where <j> = (/, A) . Without loss of generality, 
we take the kernels S 2 and S3 to be fully symmetric. As is well-known, successive derivatives of the functional W[j] 
with respect to j yield the many-body connected correlations of the field 0. 

Obviously, a first approach to handle the large- N limit of ea. (pTT|) is to use a steepest-descent method (also called 
a semiclassical or loopwise expansion in the case of usual Quantum field theory with Ti — 1/N). Thus, as described for 
instance in Zinn- Justin (1989) we may compute W[j] by expanding the action S[4>] around the saddle-point (j) s which 
yields up to order 1/N' 2 : 



W\j] = j.()) a - S((f> s ) - —TtIuGq 1 + —S 3 (x 1 ,x 2 ,x 3 )S 3 (y 1 ,y 2 ,y 3 ) 



— TrlnG^ 1 + —. h , 
2N N 2 

—G (x 1 ,y 1 )Go{x 2 ,y 2 )G (x 3 ,y 3 ) 



^G (x 1 ,x 2 )G (yi,y 2 )G (x 3 ,y 3 ) 



-0(1/N 3 ), (32) 



where we introduced the saddle-point <p s (J) and the </> s — dependent propagator Go defined by: 
<SS 



5<p(x) 



5 S 

j(x) and Go 1 (x,y;<f> a ) 



8(f>{x)8(t>{y) 



= S 2 + S 3 .<f) s , whence — -f- = S 3 , ■ JT yt = 0. (33) 



Next, we define the Legendre transform T[<p] of the functional W[j] by: 

SW ST 

T[<p]=-W\j]+j.(p with (p(x) = -—- T whence j(x) = — — . (34) 

As is well-known, the mean <f> and the two-point connected correlation G of the field <p can be obtained from T[tp] 
through the variational equations: 



(5r 

<p at the saddle-point of T: — 

dtp 



s 2 w ( 6 2 r \ 

0, and G = Ni^xt^ix^c = jj^j = [y1~J &t * = (35) 



These relations are obtained from ea. H34|) . after we notice that the physical averages are obtained when the external 
source is zero: j = 0. Note that for large N the actual two-point connected correlation {4>4>)c — G/N vanishes. This is 
natural since from ea. H31(l one indeed expects the field <f> to be "frozen" to the minimum <f) a of the action. Next, from 



the expansion (|32|l and the definition 1)34(1 we obtain the 1/N expansion of T[ip}. This yields (e.g., Zinn- Justin (1989) I: 
r[^s] = S(<p) + ^TrlnG^ 1 - S 3 (x 1 , x 2 , x 3 )S 3 (yi, y 2 , y s )G (xi, yi)G (x 2 , y 2 )G Q (x 3 , y 3 ) + 0(1/N 3 ), (36) 

where Go (a;, y; ip) is defined as in ea. H33|l but taken at the point ip. Of course, one can check in ea. l|36|l that only one 
line irreducible (1PI) diagrams 1 appear in the expansion of T[ip]. Besides, the expansion over powers of 1/N is actually 
a loopwise expansion and in ea. (|36|l we have obtained the one-loop and two-loop diagrams. Next, from ea. l|36|) and 

1 A graph is said to be "one line irreducible", or "1-particle irreducible" (1PI) in the language of Particle Physics, if it cannot 
be disconnected by cutting only one line. Otherwise, it is "1-particle reducible" . 
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ea. <|35|) we obtain the first two correlations 4> and G at the corresponding order over 1/N. Firstly, the saddle-point 
condition reads at order 1/N 2 : 



SS 



1 



5ip{x) 2N 



Tr 



6ip(x) 



AN 2 



S3S3 I G 



SGn 



5tp(x) 



.G I GoGq = at the point <p = <j>, 



(37) 



where we did not write the explicit arguments to shorten the expression. Secondly, ea. l|35|) yields for the two-point 
correlation G: 



G 1 — G 1 — S, 

where we introduced the self-energy E given at order 1/N 2 by: 



(38) 



2N 



6G 1 5G 1 

.Lrn. — z .(_ro 



2N 2 



Sip 
S3S3 Go 



Sip 



1 



2N 2 



S3 S3 G 



SG n 



S<p 



.Go 



SG n 



Sip 



Go Go Go 



SG n 



Sip 



Go G 



SG n 



S<p 



Go Gq. 



(39) 



To derive these equations we used the fact that 5 2 G n ~ 1 / SipSip — since the action S[<f>] is cubic, see eq.(j22Jl- Thus, 
once we have obtained the generating functional of the proper vertices T[tp] up to the required order over 1/N, as in 
ea. H36(l . we derive the mean <f> and the two-point correlation G from ea. (|35|) . as in eas. (|37|l - H39(l . 



4.2. Application to gravitational clustering 

The usual application of this method to a system with N fields, as in ea. (|30|l (or in usual quantum field theories), is 
to first integrate over the N fields f a since the action S n is indeed Gaussian over f a . Then, the partition function Zn 
is given by a path-integral over the one field A with a new action which contains an explicit iV-dependence. Therefore, 
we recover an expression which is similar to ea. (l31ll . to which we can apply the method described above. However, 
we cannot apply this procedure to our case since X.K S is not invertible (indeed X.K s f 2 = for all fields / such that 
/(k, w = 0, t) = 0, see eq.{7|). Therefore, integrating over the fields f a is not well suited to our problem since it would 
introduce singular terms for N = 1. Nevertheless, we can still apply the method described above to the actions given 
in eas. l|26|) - (|3U|l . Thus, in the case of the gravitational dynamics described by the action l|26|l we write: 

■ * - dlo) ■ - - d e ° - (s «' - {% R °°) • m 

where we note x = (x, p, t), the transposed matrix Rq(x, y) — Ro(y, x), and we used ea. (|23[l . Indeed, from ea. ()26jl we 
can see that N can be absorbed in the normalization of the field A and the matrix A;, hence the vanishing of the first 
two moments of A holds for any N and is recovered at any order of the expansion over 1/N. This can also be checked 
explicitly. Then, from ea. <|26|) we obtain for the matrix propagator Qo(x,y;4>): 

r.-l(„ «,7,\- ( -2X{u)K s (u;x,y) 0(y, x) - 2K s (y; x, u)J(u) \ , 
y ° {X > y ^ ) -{0(x, y )-2K s (x;y,u)7(u) -Ai(x,y) J' (4ij 

From the definition of the inverse we have: 5ijSo(x —y) = Q^.l k {x, z).Qo : kj (z, y). Substituting the expressions (|40|l and 
H41|). with A = 0, we obtain: 

0(x, z)G (z, y) = 2K s {x- z, u)J(u)G (z, y), (42) 
and: 

Ro(x, z)0(z,y) = Sd(x - y) + 2R (x, z)K s (z; y,u)f(u), 0(x, z)R (z, y) = Sd(x - y) + 2K s (x; z, u)f(u)R (z,y). (43) 

Here we used R^{x,y) = Ro(y,x) and we took the limit U — > so that Aj — > 0, see ea. (|l(J|l . Indeed, the inverse Ar 1 
does not appear in these equations so that the limit t{ — > is straightforward, which was not the case for the Feynman 
variational method l|25|l . The two eas. <|43[) are redundant but for numerical purposes it can be useful to make sure they 
are both satisfied. Thus, eas. (|42|l - H43|) yield the auxiliary propagators Go and Rq. Next, we write the actual two-point 
correlation matrix Q and the self-energy 4> as (with £(aii, X2) oc 9{t\ — t?))'- 

Q=( G R t R q) with G(x,y)=N(f(x)f(y)) c , R(x,y) = N{f(x)X(y)) c , and $ = f ° ^ ) . (44) 
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Then, convoluting from the right with Q the Schwinger-Dyson ea. ((38(l . Q 1 = Q 1 — <£>, we obtain: 

0(x, z)G(z, y) = 2K s {x; z, u)J{u)G{z, y) + £(z, z)G(z, y) + U(x, z)R T {z, y) (45) 

R(x, z)0(z, y) =S D (x-y) + 2R{x, z)K s {z; y, ujf(u) + R(x, z)E(z, y) (46) 

0(x, z)R(z, y) = 6 D (x-y) + 2K s (x; z, u)J{u)R{z, y) + £(z, z)R(z, y) (47) 

Here we have again taken the limit U — » 0. Note that eas. 1)45(1 - (147(1 are exact. To complete our system we simply need 
to add the saddle-point condition 1(35(1 and the expression of the self-energy One can easily see that at tree-order for 
the generating functional T[ip] we merely recover the linear regime ((10(1 . Therefore, we must go at least up to one- loop 
order to make progress. From ea. 1(361) this yields r = S + ^Tr In Qq . Then, the saddle-point condition 1(35(1 now 
reads: 

7^-=0: O(x,y)J(y) = ±-K B (x;y,z)G (y,z) using K s f = with /(k, w = 0, t) = 5 -^. (48) 

Here we used eq.©. The second condition ST/df is trivially satisfied. Finally, the self-energy is given at one-loop order 
by the first term in the r.h.s. in ea. ((39(l which yields: 

4 

Z(x,y) = — K s (x;x 1 ,x 2 )K s (z;y 7 z 2 )R (x 1 , z)G (x 2} z 2 ) (49) 
2 

H(x,y) = — K s (x;x 1 ,x 2 )K s (y;y 1 ,y 2 )G (x 1 ,y 1 )G (x 2 ,y 2 ) (50) 

We can check in ea. (14911 that Yj(x\,x 2 ) oc 6(t\ — t 2 ). Thus, together with the boundary conditions set by the linear 
regime (in the limit t — ► 0), eas. ((42l> - ((50l) . provide the mean /, the two-point correlation G and the response R, up 
to one-loop order within the direct steepest-descent method. We simply need to solve these equations after setting 
N = 1. One can also compute higher-order correlations through this method. Therefore, we have shown that despite 
the singular behaviour of the action S[f] in ea. 1(20(1 . which prevented the use of the Feynman variational method ((25(1 . 
we can apply the large— N expansion given by the direct steepest-descent method. Moreover, it is clear from ea. (126(1 
that all symmetries of the problem remain valid for any N and hold at any order of the expansion over 1/N. In 
particular, we recover the exact properties 1(23(1 - 1(24(1 . 



5. 1PI effective action 



Although the direct steepest-descent method described in the previous section can provide some useful insight into 
the physics at work, it has been noticed in several cases that it could also fail completely. Thus, it may exhibit secular 
terms and high-order corrections may blow up with time (e.g., Mihaila et al. (2001) \. This has led to a shift of interest 
towards other approaches to large- N expansions. We shall present in SectHJlsuch a scheme which has already been 
shown to provide a significant improvement over the direct steepest-descent method in several cases, where it was seen 
to be free of secular terms ( Mihaila et al. (2001) I or to exhibit the looked- for relaxation towards thermal equilibrium 
( |Berges (2002)||. H owever, we shall first present an alternative approach based on the "1PI effective action" T[ip] 
defined in ea. $64$ . The method described in the previous section was to firstly derive a series expansion over 1/N 
for W[j], secondly obtain through a Legendre transform an expansion for T[ip] and thirdly compute the means <j) and 
G from this expression of the 1PI effective action. An alternative approach is to firstly write the exact equations of 
motion which determine T[tp], secondly use a large- N approximation for these equations and thirdly derive (j) and G. 
That is, we directly apply the large- N expansion to T[ip) rather than W[j}. Since W and T are related through a 
non-linear Legendre transform, both procedures are not identical, although they must yield the same results up to 
the highest-order over 1/Af of the required expansion (i.e. they only differ by higher-order terms). As described in 
standard textbooks (e.g., Itzykson & Zuber (1980) I the Schwinger-Dyson equations can be obtained by noticing that 
the integral of a derivative vanishes: 



54>{x) 



,N\j.4,-S{4>)] 



whence 



J'0) 



SS 



S<j)(x) \N Sj 



1 5 



■m = o. 



For the cubic action S[<j>] defined in ea. 1311) this yields: 

5W SW 



SW 1 

j(x) - Sl(X) - S^X^Xi).— r - -S 3 (X,X!,X 2 ) 

dj[xx) 2 



1 



S 2 W 



5j{xi)8j{x 2 ) N 6j(xi)5j(x 2 ) 



(51) 



(52) 
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By taking successive derivatives over j at the point j = of eq. i|52|) we would obtain a hierarchy of equations between 
the many-body correlations G q . In this way we would recover the standard BBGKY hierarchy (Peebles (1980) I. 



However, here we introduce the 1PI effective action defined in ea. ((34(l and we write ea. ((52(l in terms of T[ip\. This 
reads: 



5 JL _ 8 A _ _L S G-0 with G - — - ( \ 1 (53) 
6(p Sip 2N 8j6j \6(pS(p J 

As seen from ea. (|35H G/N is actually the two-point correlation of the field <fi at the point ip = </>. The equation 1(53(1 
is exact. We only used the fact that the action S[c/>] is cubic as given in ea. l|3*l"|) . From eq.JJU we note that at the 
saddle-point <fi of T we have the relation: 

5 S 1 

— + — S 3 .G = at the saddle-point ip = 6 and with G = N{4>4>) c . (54) 
dip 2N 

The ea. 1(53(1 is an intricate non-linear equation for the functional T[ip] which we do not know how to solve. However, 
by taking successive derivatives of ea. ((53() with respect to tp we can obtain a useful hierarchy of equations. A first 
derivative of ea. l{53")) yields: 

G- = G --E with ^( W ) = ^^ and E(x,y) = -±S 3 (x, ., .)|^. (55) 

The propagator Go and the self-energy E are defined as in ea. l|33|l and ea. (|38|l . However, the last expression for E in 
ca. l|55|l is exact while ea. (13911 is only its expansion up to order 1/N 2 . Introducing the three-point vertex T 3 we can 
also write: 

1 (5 3 r <5G _1 

E(x,y) = — S 3 (x,ui,U2)G(ui,u 3 )T s (y, u 3 , u 4 )G(u 4 , u 2 ) with r 3 = $ $ § = ^ ■ (56) 

Next, taking a second derivative of ea. l(5*5|) we get: 

T 3 -S 3 -^-S 3 .G.T 3 .G.T 3 .G+^-S 3 .G.T A .G = with r 4 = /T . (57) 
N 2N dipoipoipdip 

Obviously, by taking successive derivatives of the equation of motion 1(53(1 we obtain an infinite hierarchy of equations 
which link the successive proper correlations T q . However, because of the factor 1/N in ea. (|53fl higher-order vertices 
T q contribute to the self-energy E at a higher order over 1/N. Therefore, we can close this hierarchy by looking for an 
expression of E up to a given order over 1/N. Moreover, we note that since the action S[cj>] is cubic there is no "bare 
contribution" S q to the q— point vertex T q for q > 4 so that T q ~ 1/N for q > 4. Note that this is consistent with 
ea. 1)36(1 . Hence we obtain at order 1/N for T 3 : 

T 3 (x, y, z) = S 3 (x, y, z) + j^S 3 (x, u 1 ,u 2 )G(u 1 ,u 3 )T 3 (y, u 3 , u 4 )G(u A , u 5 )T 3 (z, u 5 , u 6 )G(u 6 ,u 2 ) + 0(1/N 2 ). (58) 

Thus, we obtain a closed system and we can solve for <fi, G and T 3 through eas. l|54l) - (|58|l . Of course, the equations 
derived in this section agree with those obtained in Sect. 0]up to order 1/N 2 since they were both derived through 
1/N expansions up to this order. However, they differ by higher-order terms. This is obvious from the fact that E 
now depends on the exact propagator G so that solving exactly the equation ((55(1 , with E obtained from eq. (15811 , will 
generate contributions to G at all orders over 1/N. In other words, by going from the method presented in Sect0] 
to the present method we have made some infinite partial resummations. This also means that we must solve these 
equations exactly. Indeed, there is no point to solve them as a perturbative series over 1 /N since this would give back 
the results obtained from the direct steepest-descent method described in Sect 01 

Since both this 1PI effective action method and the former steepest-descent approach are derived up to the same 
order over 1/N one might naively expect their accuracy to be of similar order. Moreover, although the 1PI effective 
action scheme includes partial resummations with respect to the other method, it is not obvious a priori that this 
should improve the final results. Indeed, we still miss some of the higher-order terms and we can easily imagine peculiar 
cases where a resummation of some badly chosen diagrams would actually dramatically worsen the results (see also 
Sect lG.4l below'). However, these new equations of motion for / and G were seen in several instances to yield a significant 



improvement over the steepest-descent scheme (e.g., Mihaila et al. (2001) I. This can be understood in a more generic 



context from the fact that the self-energy E now depends on the actual correlation G so that the Schwinger-Dyson 
equation G _1 = Gq 1 — E is now an implicit equation for G. Since implicit schemes are usually more stable than 
explicit ones and satisfy physical constraints it is not surprising that the 1PI effective action approach should perform 
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better (a familiar example is provided by the numerical integration of ordinary differential equations). This can also 
be analyzed in the light of virial scalings for the case of the gravitational dynamics, see the discussion in SectlHJ 

Finally, we note that these equations were also used in Mihaila et al. (2001) to study the quantum roll. Going to 
order 1/N for S, they called this approximation the "bare vertex approximation" since at this order we simply have 
T 3 = S3. Here we have seen that one can actually go to any order over 1/N by including all vertices T q up to the 
required order. The last vertex is indeed given by its "bare contribution" S q , which in our case yields T q — if we close 
at a level q > 4. Thus, we see that we have actually closed the infinite hierarchy of equations for the correlations T q by 
neglecting the high-order proper correlations T q beyond some finite order. This is somewhat similar to the standard 
closure of the BBGKY hierarchy. However, the truncation of the proper correlations T q is much more powerful since it 
should still be able to handle the non-linear regime. This is not surprising since, as is well-known, as soon as 1^3 7^ one 
obtains a non-zero estimate for all many-body correlations G q . Therefore, working with T[<p] allows us to investigate 
powerful closure schemes which may still handle the strongly non-Gaussian regime. Note also that this 1PI effective 
action method provides an estimate for all many-body correlations G q (since from T we can derive W) although in 
this paper we restrict ourselves to / and G. 



6. 2PI effective action 

6.1. General formalism 

In the 1PI effective action approach developed in Sect|S]one needs to study higher-order proper correlations T q as 
one pushes the method to higher orders over 1/N. This is not very convenient for numerical purposes since it implies 
that one needs to solve implicit non-linear equations for quantities which depend on an increasingly large number 
of variables. This is especially critical in our case where the coordinate x = (x, p, t) is actually seven-dimensional so 
that one very quickly goes beyond computer power. Therefore, we introduce in this section another large— N method 
which is likely to be more powerful to study large-scale structure formation. It only involves the one and two-point 
correlations <j> and G, whatever the order over 1/N, but it again yields an implicit equation for G (which is a key 
feature to handle the non-linear regime). This can be done by introducing higher-order Legendre transforms (e.g., 



De Dominicis (1962) |, which are a direct generalization of the 1PI effective action T used in Sect|S| Thus, following 



Cornwall et al. (1974) we define the double generating functionals Z[j,M] and W[j, M] by: 

Z[j, M] = e w ^ M ' = J [64>] e ^W+HM.0-s(0)] ; (59) 
for any test-field j(x) and symmetric matrix M(x\,X2). Next, we define the double Legendre transform T[ip, G] as: 

(60) 



1 1 SW sw 1 

r[<p,G} = -W\j,M]+j.<p + ~<p.M.v+—Tr[G.M], with — = <p and m = ^ 



1 „ 

(pip + —G 



N 



This is the generalization of eas. i|31|) . i|34|) . In particular, as in ea. l|34|l we also have: 



$r $r 1 $r 

— = j + M.tp and — = — M, while F[<p] = TUp, G] at the saddle-point G such that — = 0. (61) 
S<p ' SG 2N ~ SG 

The last relation in ea. l|61|l gives the 1PI effective action T[ip] used in Sect|S]in terms of the new double Legendre 
transform T[p, G}. Besides, as in ea. ()35|) the mean <j> and the two-point correlation {4>4>)c are obtained from the 
variational principle: 

1 (5r (5r 

<p — tp and (</>0} c = — G at the saddle-point (ip, G) such that — = 0, — = 0. (62) 
N Sp SG 

Moreover, while T[p] defined in Sect^was the generating functional of 1PI diagrams with "bare propagator" S^ 1 
the double Legendre transform T[<p, G] is the generating functional in ip of two-line irreducible (2PI) diagrams 2 with 
propagator G (e.g., |De Dominicis fc Martin (1964)[ Vasil'ev fc Kazanskii (1972)| see also Apd Ia")) . Then, following 



Cornwall et al. (1974)| it is convenient to write T\ip, G] in the form 



I>, G] = S[<p] + ^Tr [G^.G - l] - ^Tr In (A" 1 ^) - if fc>, G] with A -1 = S 2 , G^ 1 = (63) 

2 A graph is said to be "two line irreducible", or "2-particle irreducible" (2PI), if it cannot be disconnected by cutting only 
two lines. Otherwise, it is "2-particle reducible". 
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Note that the propagator Go is defined as in eas. (|33[) and (|55|l while A is the "bare propagator" of the action S [</>]. 



Thus, ea. H63|) defines the new generating functional f. Then, as shown in Cornwall et al. (1974) the reduced 2PI 



effective action T[<p, G] is actually the sum of 2PI vacuum graphs of a theory defined by the action S[(f>] given by: 

S[(f>] = ^-(fj.G^ 1 .cf) + NSi n t[4>; <f] where S- m t[(f)\ <p] = cubic and higher-order terms over <p of S\ip + (j>\. (64) 

As seen in ea. (|64|l . in general the action S[<j>] depends parametrically on the field ip, whence the dependence on ip of 
the functional T. However, in the case of a cubic action as l|31|l . we see that the interaction part S- m t[4>; <p] introduced 
in ea. Q64[l does not depend on ip. Therefore, the functional T does not depend on ip either and we have: 

for a cubic action S: f [G] = 2PI vacuum diagrams of the action S[<p] — —<j).G^ 1 .<j) + — S 3 (f> 3 . (65) 



We give a simple proof of the result l|65|l in App^J following the analysis used in Vasil'ev fc Kazanskii (1972)| to show 



that the full effective action T[ip, G] only contains 2PI diagrams. Indeed, the case of a cubic action S[<p\ where the 
reduced 2PI effective action T only depends on G is much simpler than the general case. Since this is actually the 
relevant case for gravitational dynamics and the formalism associated with the double Legendre transform T[ip, G] is 
not so widely used, Add lAl makes this article self-contained. In particular, the analysis presented in Add 1X1 should be 
sufficient for readers who are only interested in the application of this method to the problem we investigate here: the 
formation of large-scale structures through gravitational instability. The result (|65|l provides a loopwise expansion of 
r[G] over powers of 1/N. From this expression for T, we can derive the mean cj) and the correlation G through the 
variational equations l|62|). Besides, from F we can also obtain W and all correlations at any order. Of course, we can 
see that the saddle-point condition 162|) for G yields a non-linear implicit equation for G (since T is a functional of G) 
as we wished. More precisely, let us write ea. (|6*2")) using ea. (l6*3"j) . This reads: 



5T _ 5S 1 
~fcp~ ~ ° 1 + 2iV Tr 



Sip 



SS 1 

whence — + — S 3 .G = 0. (66) 
dip 2N 



Thus, we recover the exact equation Q54JI ■ Indeed, so far we have made no approximation. We only used the fact that 
the action S[<j>] is cubic. Next, the saddle-point condition for G writes: 



§ = 0: G-=G-- S with S = 2§, 
This equation gives the exact expression of the self-energy S in terms of T. Note that ea. l|67|) clearly implies that T 



^ = 0: G- l =G^-V with S = 2^. ((i7) 



must be two-line irreducible (2PI) since, as is well known, only 1PI graphs contribute to E, see also Berges (2002) 



6.2. Comparison with the 1PI effective action approach and the BBGKY hierarchy 

In order to compare this method to the 1PI effective action scheme described in SectEl let us compute the two-loop 
and three- loop contributions to T. This yields: 

^ = j^^SsG 3 + -^j^S3S 3 S 3 S 3 G 6 , (68) 
which gives for the self-energy E: 

%(x,y) = -^S 3 (x,x 2 ,x 3 )S 3 {y,y 2 ,y 3 )G(x 2 ,y2)G(x 3 ,y 3 ) 

+ 2jj2 S a( x > x %i x 3)S 3 (y, j/2, y 3 )S 3 (z 1 ,z 2 , z 3 )S 3 {u u u 2l u 3 )G(x 2 , z 1 )G(x 3 ,u 1 )G(y 2 , z 2 )G{y 3l u 2 )G{z 3 , u 3 ). (69) 

Then, we see that both methods are actually identical when we close the hierarchy in Sect|S]at order 1/iV and we only 
keep the two- loop term in ea. (|68|l . However, when we go to higher orders (here 1/A^ 2 or three-loop) these two methods 
are no longer identical, but of course they still agree up to the order over 1/N up to which they were derived (i.e. 
they differ by higher-order terms). This is not surprising. Indeed, as described in Sect|SJ in the previous 1PI effective 
action method as we go to higher orders over 1/N we close the hierarchy of equations for the vertices T q at higher 
levels. Hence we include an increasingly large number of exact equations of motion for these correlations T q (hence 
for the connected correlations G q ). This can be a nice feature, but as we noticed above it is not very convenient for 
numerical purposes since it leads to non-linear implicit equations which involve higher-order correlations. By contrast, 
in the approach discussed here the basic unknowns are the field ip and the correlation G, whatever the order over 1/N. 
Indeed, although high-order diagrams for T can lead to lengthy expressions for the self-energy E they only involve the 
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explicit vertex S3 and the two-point correlation G. Of course, the drawback is that the only exact equation of motion 
of the previous hierarchy which remains valid within this approach is the first relation (|66|l . 

This also clearly shows how we could extend this method to include exactly some higher-order equations 
of the hierarchy. We simply need to work with higher-order Legendre transforms T[tp, G, G3, .., G q ] (see also 



Calzetta & Hu (1994) I. Then, we would recover all exact BBGKY equations between the first few correlations 



6, G, G3, .., G q , up to order q. Of course, as in Sect[51 the advantage of the method described here with respect 
to a simple truncation of the BBGKY hierarchy itself, is that higher-order correlations G 9 +i,.., are not taken to 
be zero. Rather, any expression for the functional V (or T) defines a precise closure where higher-order correlations 
(G 9 +i, ..) are expressed in terms of the lower-order correlations (<fi, G, .., G q ) which are explicitly taken into account 
in the Legendre transform T. Thus, by working with the double Legendre transform T[ip, G] introduced in ea. (|60() . 
we keep the exact BBGKY relation between cf> and G and we approximate higher-order correlations G3, .., by some 
functionals of <f> and G. Of course, this also means that these higher-order correlations G q will no longer satisfy the 
exact BBGKY hierarchy. This is why the 2PI effective action method described here is not identical to the scheme 
presented in Sect[S] 

However, once we have an approximate 2PI effective action T[ip, G], which also defines the usual 1PI effective action 
r[ip) through ea. (|61|l . then all many-body correlations G q will satisfy the new BBGKY hierarchy associated with this 
new physical system (which we define by its 2PI effective action). Thus, the method described in this section is rather 
elegant as it keeps the variational structure l|62|) and the Legendre relation l|f)()|l . Besides, in the context of Quantum 
field theories it can be shown that in some specific cases (Baym (1962) I some very important properties can be proved 
from the mere fact that the equations of motion can be derived through a variational principle as in eq. I|67|) . In our 
case, these consequences have no practical interest but it might happen that another similar effect comes into play 
(note indeed that the fact that the self-energy is the derivative of some functional 2r[G] is a non-trivial constraint, 
especially when the field 4> is actually a composite field as in our case). 



6.3. Equations of motion for gravitational clustering 




Fig. 1. The two- loop 2PI vacuum diagrams of the reduced effective action T. Each solid line stands for a propagator 
G, each dashed line for the propagator D and each full-dashed line for the mixed propagator R (cross-correlation). 
The indices /, A denote the /—field and A— field sides of the propagators G, R or D. These indices /, A also correspond 
to the legs of the three-leg vertex K s (big dots). 



We now apply to gravitational clustering, described by ea. l|26() . this 2PI effective action method. Thus, ea. (|66|l 
reads: 

^r=0: 0(x,y)J(y) = ±-K s (x;y,z)G(y,z), (70) 

while the second condition ST/5f is again trivially satisfied. Note the difference between eq.{7UJ) and ea. (|48|l obtained 
in the direct steepest-descent method. As explained above, ea- lfTOjl is actually exact. In fact, for N — 1 it can be directly 
obtained by taking the average of the Vlasov equation (JSJ (with ti — > 0). Next, the Schwinger-Dyson equation l|67|l 
yields again eas. (|45|l -H47 |l as in Sect0] However, the self-energy $ is different and it is now derived from the functional 
T, which is obtained through a loopwise expansion of the 2PI vacuum diagrams defined by the action S. The latter is 
given by: S- m t[4>] = —^-K s f 2 and we now have three propagators: G,R and D — N(\X) C . Since D eventually vanishes 
it was not needed within the direct steepest-descent scheme presented in SectQJ However, within the 2PI effective 
action method we must keep all terms in the intermediate steps. Then, the two-loop contribution to T is given by the 
two graphs displayed in Fig^ which yields: 

f 2 ioo ps = ^((-\.K s ff)(-X.K s ff)) 2Pl 

= j^K s (x;x 1 ,X2)K s (y;y 1 ,y 2 ) [D(x,y)G(x 1 ,y 1 )G(x 2 ,y 2 ) + 2R T (x, y 1 )G(x 2 , y 2 )R{xi, y)] . (71) 
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Next, from eq.J7TJ and ea. l|67|) . which now reads $ = (2ST /6Qij), we obtain the two- loop contribution to the self-energy. 
This yields: 

4 

^2\oo ps (x,y) = —K s (x;x 1 ,X2)K s (z;y,Z2)R(xi,z)G(x2,y2) (72) 
2 

n 2 ioo P s(^, y) = —K s (x;x 1 ,X2)K s (y;yi,y2)G(xx,yi)G(x2,y2) (73) 

Moreover, we can check that D and the (1,1) component of $ vanish. Note again the difference between eas. <|72l) - H73l) 
and eas. H49l) - (|50|l . By going to higher orders over 1/N (higher-loop diagrams) we would add higher powers over f,G 
and R to the self-energy. We obtain a finite number of terms because to each order over 1/N only corresponds a 
finite number of diagrams in the expansion of the reduced 2PI effective action T. This nice feature arises because we 
kept the auxiliary field A throughout the calculation. Since the actions in eas. (|26[l - (|30|> are quadratic over A we could 
have performed the Gaussian integration over A. Then, we could apply the same formalism. However, the drawback 
of this approach is that to each order over 1 /N would now correspond an infinite number of diagrams which we would 
need to resum (see Aarts et al. (2002)] ) . Thus, keeping the auxiliary field A allows us to perform these resummations 



in an automatic manner through the matrix R. Moreover, as shown in Sect l3.ll the field A actually generates the 
response-functions of the system so that R also has a physical interpretation as the mean first-order response function. 



6.4. Infrared divergences 

Although we applied the 2PI effective action method to the action S[f, X] given in ea. l|26[) . we could also apply it 
to the (N + 1)— field action H30|) . looking for a symmetric solution. The analysis is somewhat more intricate because 
we now have four different propagators G, G, R and D (G corresponds to {f a fb)c with a ^ b while G corresponds 
to (f a fa)c)- Moreover, these propagators now exhibit different scalings over N. Thus, there is no longer a one-to-one 
correspondence between the loop-order of a contribution to the reduced effective action T and the order over 1 /N of 
the required correlation G. One needs to check explicitly by inspection of the equations of motion the order over 1/N 
required for each self-energy term to obtain the correlations G, G, R and D up to the looked-for order over 1/N. In 
particular, within this framework the diagram to the right of Fig^ yields a subleading contribution as compared to 
the left diagram so that it would be disregarded. This implies that the frameworks (|26J) and (1.>0I) are not equivalent 
for JV ^ 1. 

However, before we can apply any of these procedures (i.e. working with (|26|l or (J20J0 we must pay attention to 
the long- wavelength divergences which we discussed in Sect 12. 31 We recalled that within the framework of standard 
perturbation theory these infrared divergences cancel out, at least at leading order (e.g., |Jain &: Bertschinger (19 96)). 
Obviously, we need to preserve this property in our approximation scheme. Otherwise, our results would be grossly 
inaccurate and useless (even though realistic power-spectra have n > — 1 on scales x ii 10 Mpc, this effect would 
prevent accurate results at smaller scales which are precisely those of interest). Therefore, we must make sure that we 
preserve these delicate cancellations between different diagrams. Thus, once we have selected all diagrams which are 
needed to reach the order over 1 /N which we wish to obtain, we must add the sub-leading diagrams (if any) which 
contribute to higher orders over 1/N but provide the right "counter-terms" (at N = 1) to the infrared divergences of 
the graphs we have already included. 

The identification of the right diagrams to group together in order to cure these infrared problems is not obvious 
in the expansion over 1/N if we work with the (N + 1)— field action l|3U|) . Indeed, it is a completely unrelated point. 
Fortunately, this question is easily solved by working with the "semiclassical expansion" (|26l) . Indeed, as seen for 
instance from the second expression in ea. (|26|) we see that the factor N can be absorbed within the normalization of 
the correlation A;. Therefore, the expansion over 1/N is actually an expansion over the amplitude of Aj, hence over 
powers of the linear power-spectrum PLo(k), see Sect l2.2l As a consequence, it shows the same cancellations of the 
infrared divergences as the usual perturbative approach for any N and at any order over 1/N. On the other hand, 
since the (N + 1)— field approach does not affect the same powers of 1/N to the various diagrams it does not exhibit 
this infrared cancellation (except for N = 1 if we sum all diagrams). Therefore, it is not well-suited to gravitational 
clustering studies 3 , unless we use the vertex given by the modified Vlasov equation Ijl4|) which automatically damps 
long-wavelength contributions. However, as noticed in Sect l2.3l this implies the loss of homogeneity. Hence, we can 
conclude that to investigate gravitational clustering we had better use the "semiclassical" expression l|26[) rather than 
the (N + 1) -field generalization (3DJ). 

3 It is interesting to note that in a recent paper, studying for a 1 + 1 dimensional (j> 4 field theory several large- N schemes like 
those described in this article, [Cooper et al. (2003)| found that the accuracy of the numerical calculations as compared with the 
exact result at N = 1 was better when one keeps both diagrams in Fig0 We have seen that this is also the case here because 
of the infrared problems. One may wonder whether this is a general feature. 
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7. Perturbative expansion 

Finally, we now solve eas. H42(l - (|50[l as a perturbative series over the linear power-spectrum Plq- This actually boils 
down to the standard perturbative approach and exhibits the same ultraviolet divergences. Hence it cannot handle the 
non- linear regime but our purpose here is to check that the infrared divergences discussed in Sect 16.41 indeed disappear 
when we compute the equal-time correlation of the density field. In this section, we consider a critical-density universe 
and we normalize the time t and the scale factor a to unity at z = 0. We also define a new time coordinate r by: 

t =- ln(</t ) = - lni, hence t = e 3r , a = e 2r and p = — !— . (74) 
3 3 QttQ 

Then, the linear operator O and the vertex K write in Fourier space: 

d d k w 

O = — - 3e~ r k. — , K(x; x u x 2 ) = (2^) 3 2e T '<5 D (n - t)8 d (t 2 - t)5 d {V. x + k 2 - k)5 D (wi)fe(w 2 - w)-^p, (75) 

while from Sect 12. 2| the two-point correlation = (??i(ki,Wi, ri)rji(k2, W2, T2)) c of the linear growing mode reads: 

(76) 



A i (k 1 ,w 1 ,r 1 ;k 2 ,w 2 ,r 2 ) = fe(ki +k 2 ) P L0 (*i) e 2Tl+2T2 

(Z7T J u 



3 fc 



2 



2 ,, k 2 .W2 
1 + -e T2 — — - 



3 k 



2 



At lowest order, the mean / and the propagators Go and G obtained from eas.l|42 |) -l|5U |) are: 

J i0 \x) = ^ 6 D (k), G^\x 1 ,x 2 ) = G^(x 1 ,x 2 ) = A L (x u x 2 ). (77) 

This simply corresponds to the linear regime. On the other hand, from eas. (|43|) - (|47|) we obtain for the response function 
R = (/(ki,wi,Ti)A(k2,W2,T 2 )) c : 

R { °\x u x 2 ) = R^(x u x 2 ) = ^-g fe(ki + k 2 ) 6(n -r 2 )| ^(w! + w 2 + 3k!(e- T2 - e^ 1 )) 



- p dr Mw 2 + 3k!(e-^ - O) (~e T ^Ji [2 e 3 ^ 3 - + 3 e -^+ 2T ] + | [ ( 
J T2 t O Kj 5 



We can check that it is indeed causal and it satisfies both eas. (|43|l . Moreover, it also obeys both exact constraints 
124(1 . In order to check the second condition 124JI one must pay attention to the large-scale cutoff of the gravitational 
interaction which multiplies the factor 6/5 in the last term by k\j(k\ + fc 2 ) (see eq.JSJl). At next order, we obtain from 
ca. 1481) for the mean /: 

7 (1) (x) = K s (x;y,z)G^(y,z) = 6 D (k) 1 J dk' P L0 (k>) (^^j , (79) 

where K s is the integrated vertex Q1H)). Note that the integral over k' diverges at long wavelengths for n < —1. This 
corresponds to the infrared divergences discussed in Sect l2.3l and Sect 16.41 Next, in order to compute the next-to- leading 
contribution G^ to the correlation G we first obtain the self-energy at lowest order from eas. l49|) - (l5U(l . Substituting 
these expressions into ea. <|45|) we can derive G^ 1 '. Here we are only interested in the density fluctuations hence we 
define: 

G(ki,0,n;ka,0,7n)=<y £ ,(ki+ka)G(fci;Ti s 7a) whence P(fc, t) = (2^) 6 G(fc; t, t), (80) 

where P(fc,r) is the non-linear power-spectrum. Thus, G{k]T\ 1 T 2 ) yields the two-point correlation of the density field 
(in spatial Fourier space) at times n and t 2 (indeed w = corresponds to the integration over the momentum p). 
After a tedious calculation, we obtain from ea.Q45[l: 

GW(fc;r 1 ,r 2 ) = 9(n - r)6(l - e T ^)G^ (fc; r, r 2 ) + -l-P LQ (k)P L0 (k')e^+ 2 ^ ^-F^(k', -k', k) 

P LO (k-k')PLo(fc')e 4Tl+4T 4 i? 2 (s) (k',k-k') 2 , (81) 



is) is) 

where the symmetric kernels Fk and Fi are those which appear in the standard perturbative expansion derived 



2tt) 6 uv ' x ' 5 
.vuac uiic ojiUMVMnC kernels F 2 an d F \ are those wm. 

from the hydrodynamical approach (e.g., pain fc Bertschinger (1996*} I. In particular, we have: 

(82) 



pWflr' w irM- 5 i l k, -( k - k | I k'.(k-k') 2 [k'.(k-k')] 2 
2 1 ' j ~7 + 2 fc'2 + 2 Ik-k'P + 7fc' 2 |k-k'P 
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and: 

(«) rv u/ - 10fc/4fc6 + 10fe 4 fc /6 ~ 21fc 6 (k.k') 2 - 44fc 4 fc /2 (k.k') 2 - 59fc 2 fc /4 (k.k') 2 + 76fc 2 (k.k') 4 + 28fc' 2 (k.k') 4 , 
F 3 (k,-k,k)- 126fc' 4 fc 2 |k-k'| 2 |k + k'| 2 (83) 

Next, the symmetric solution of the linear equation (|81|l is simply: 

G«(fc;r 1 ,r 2 ) = ^ P LO (fc)P L o(fc')3F 3 (s) (k', -k', k) + ——P L0 (k- k')P i0 (fc')2if } (k', k- k') 2 (84) 

Therefore, we obtain for the non- linear power-spectrum up to order -P 2 Q : 

P(k,r) =e 4T P L0 (k)+e ST \p L0 (k)P L0 (k , )6F^\k\-k\k)+P L0 (k-k')P L0 (k , )2F^\k\k-k') 2 ] + O(P 3 L0 ) (85) 



We can see that we recover the results of the standard perturbative approach (e.g., Goroff et al. (1986^| 



Jain & Bertschingcr (1996)). In particular, the integration over k' in ea. H85|) no longer diverges at long wavelengths 



Therefore, the system l|42|) - (|50|) is well-behaved with respect to these infrared divergences and yields meaningful re- 
sults. However, note that we must still have a large-scale cutoff, or n > — 1 on large scales, so that the intermediate 
quantities / and G(x\,X2) are finite. Besides, we explicitly see in eg. (18411 that the cancellation of the infrared diver- 
gences only holds for the equal-time correlation of the density field. This is natural, since by looking at the density 
field at two different times t\ ^ r 2 we are obviously sensitive to the mean value of the velocity field which translates 
the small-scale density fluctuations over some non-zero distance during this finite time-interval. 

Next, we can easily see that the same analysis can be performed at next-to-leading order for the 2PI effective 
action equations derived in SectEl Therefore, both sets of equations are well-behaved with respect to the infrared 
divergences. However, note that within the direct steepest-descent approach although the actual non-linear correlation 
G shows the correct infrared cancellations this is not the case for the intermediate propagator Go, but this is not a 
problem since Go has no physical meaning. 

Finally, we can note that eq. (|84J) confirms that solving the Vlasov equation through a perturbative expansion over 
powers of the linear power-spectrum PLoik) recovers the results obtained from the hydrodynamical approach. This 



is consistent with the analysis presented in Valageas (2001) This also clearly shows that we must solve exactly the 



equations l|42|l - H5U|) or l|7U|) - (|73|l in order to obtain interesting results. As compared with the perturbative expansion 
over powers of PLo(k), this automatically performs some partial infinite resummations and it is expected to be free of 
the ultraviolet divergences encountered in the former method. Within the 2PI effective action approach, the advantage 
of computing exactly the mean / and the correlations G and R, once we have derived T[Q] up to some finite order, is 
that it preserves the variational principle (|62[) . 



8. Conclusion and discussion 

In this article we have developed a new approach to investigate the dynamics of gravitational clustering in an expanding 
universe. The main results of this paper are the following. 

- We discussed the long- wavelength divergences which may appear for some initial conditions fSect l2*3|) and we 
showed how they could be treated by a simple modification of the equations of motion. However, this leads to a loss 
of homogeneity; hence it is more convenient to keep a large scale cutoff and to make sure that any method used in 
this case shows the same infrared cancellations as the standard hydrodynamical perturbative expansion. 

- We derived the actions S[f] and S[f, A] which allow one to obtain the statistical properties of the phase-space 
distribution function /(x, p, t) from a path-integral formalism fSect !3.l|l . We noticed that the auxiliary field A generates 
the response functions of the system. 

- We briefly discussed several approximation schemes which may be applied to this functional formalism (Sect[!OJ). 
We recalled that a perturbative expansion over the non-Gaussian terms only recovers the standard hydrodynamical 
results. Then, we showed that the Feynman variational method cannot be applied to our system because of the 
singularities of the action S[f]. Next, we presented a "semiclassicai" approach and an Af— field generalization which 
may serve as a basis for large— A^ expansions. 

- We derived the equations of motion up to next-to-leading order within a direct steepest-descent method for 
gravitational clustering (SectQJ. We showed that this is possible despite the singular behaviour of the action S[f]. 

- We presented another scheme (based on the 1PI effective action) which is expected to be more powerful than the 
previous steepest-descent method by analogy with other physical problems (Sect|SJ). In particular, we compared both 
approaches (with some emphasis on explicit versus implicit equations) as well as the BBGKY hierarchy. 

- We eventually derived for the gravitational dynamics the equations at next-to-leading order given by a 2PI 
effective action approach (SectEJ- This is the most elegant of these three methods, which we compared in Sect 16. 21 
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- We showed that the infrared divergences must be taken into account and select the "semiclassical" approach over 
the N— field generalization fSect !6.4) l. Indeed, only the former scheme shows the right infrared cancellations. In fact, 
the "semiclassical" method agrees with the standard perturbative results up to the order to which it is derived (and 
it includes in addition a partial resummation of higher-order terms). 

- We checked that a perturbative solution of the equations obtained from a large- TV approach up to next-to-leading 
order indeed recovers the perturbative hydrodynamical results at this same order ( Sect 01 . We noticed that the large- 
scale regularization of the gravitational interaction must be kept in mind, and that the infrared divergences which 
occur for some initial conditions only cancel for the equal-time density correlation. 

Finally, we note that the direct steepest-descent method and the effective action approaches can be compared in 
the light of the virial scalings which are a key feature of gravitational systems. Thus, an isolated gravitational system 
obeys at equilibrium the property 2T+V = 0, where T (resp. V) is its kinetic (resp. potential) energy. In a cosmological 
context, halos are not isolated nor at perfect equilibrium but they still obey a "virial scaling" T ~ \V\. Thus, after an 
overdensity has collapsed, the physical velocity scale v obeys v 2 ~ QM/R ~ QpR 2 within an object of mass M, size 
R and density p. This "equilibrium" (even though not perfect) between the kinetic and potential energies plays a key 
role in the non- linear regime of large-scale structure formation where such virialized objects build up (e.g., clusters 
or even filaments). Then, one can see that the implicit Schwinger- Dyson equations (|55|l or i|67|) for G ensure that the 
"virial scaling" will be satisfied, provided one can attach a unique time-scale and velocity-scale to a given length-scale, 
since the self-energy explicitly depends on the non-linear two-point correlation G (of course this "virial scaling" is also 
apparent in the exact BBGKY hierarchy). Thus, at next-to- leading order, the Schwinger-Dyson eas. (|45|I - I47|I with the 
self-energy (|72[1 - (|73[1 yield through a simple dimensional analysis: 

OG ~ K S K S RGG whence 1 ~ (*fcfc 3 ^r) f, f(k,w k ) ~ t^k^w^ 1 and p(x) ~ k/(t k w k ) ~ t^ 2 , (86) 

where we assumed that the time-scale t k and the velocity-scale 1 /w k are associated with the length-scale x = 1/k 
(whence t% ~ Wk/k). The result p ~ is the "virial scaling", where p is the real space density and t k would also be 
called the local dynamical time. By contrast, within the direct steepest-descent scheme it is clear that in the regime 
where G becomes very different from Go the explicit equations (|49J) - H5()|I imply an increasingly large deviation from 
the "virial scaling". This means that the direct steepest-descent method may not be able to handle the non-linear 
regime since it generically misses the "virial scaling" . On the other hand, the effective action schemes automatically 
recover this behaviour (which is clearly a very important non-perturbative property of the system) because of their 
implicit form. Note however that this analysis does not mean that the "virial scaling" holds for any scale k. Indeed, 
as seen in (|8fif) in order to derive this property one needs to assume that the physics at scale k is governed by this 
local scale and some associated time-scale t k and velocity-scale l/w k . However, this is not necessarily the case because 
in the context of large-scale structure formation in an expanding universe there are always a characteristic time tn 
(the Hubble time) and a characteristic scale Rq (the scale which is just turning non-linear) so that the dimensional 
analysis above is not sufficient (e.g., one can introduce any dependence on the dimcnsionlcss factor kRo). However, 
we can expect the "virial scaling" to hold at scale Rq where t% ~ tn- 

Of course, the actual test of the interest of the methods developed in this paper will come from a numerical 
computation of the equations of motion derived here. This is in itself a non-trivial task, since the basic object is 
the phase-space distribution /(x, p,t) which depends on seven coordinates. We plan to investigate such a numerical 
calculation in future works. On the other hand, one may apply further approximations to the equations of motion 
derived in this article at next-to-leading order over l/N. Finally, we note that the path-integral framework developed 
in this article may also serve as a basis to other tools borrowed from field theory than the "large— TV methods" we 
focused on. 



Appendix A: Reduced 2PI effective action f 

In this appendix we derive eas. (|6^|i - (|65|l for a cubic action. That is, we show that the reduced effective action T defined 
in eq.ljnSJ) is indeed given by the 2PI vacuum diagrams of the action S[4>] defined in ea. (|65fl . The case of a cubic action 
is much simpler than the general case since T[tp, G] no longer depends on tp. Thus, we present below a diagrammatic 
proof of eas. l|rjSl) - (|631) which is much simpler than the method described in Jackiw (1974)|Cornwall et al. (1974~ We 
actually follow the analysis used in Vasil'ev & Kazanskii (1972) to show that the full effective action T[ip, G] only 
contains 2PI diagrams. Indeed, for a cubic action this method also allows one to show that the reduced effective action 
r[G] is given by a sum of 2PI vacuum diagrams. To simplify the expressions, in this appendix we take N = 1 so that 
we write: 



Z\j, M] = e w ' J ' M ' =M [dcj>(x)} ei-<H-h*-M-*-sW with S y,] = Sl . 



Va- 1 . 



(A.l) 
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The action S[(/>] is identical to eq.JjHJ with non-local terms and we defined the "bare propagator" A by A -1 = 52- 
Besides, we normalize the generating functional Z[j,M] and Z[j] (the latter being defined as in ea.ipnjl. that is 
without the symmetric external matrix M) by the factor M defined by: 

Af(A) J [&((>} e -**- A " 1 -* = 1, hence Af(A) = (DetA^ 1 ) 1 / 2 . (A.2) 

Thus, if Si = and S3 = we have Z(j = 0) = Z(j = 0, M = 0) = 1. On the other hand, both functionals are related 
by: 

Z(A;j,M) = j^-Z(A-J) with AT^A^-M, (A.3) 

where we wrote explicitly the bare propagator corresponding to each generating functional Z and we introduced the 
new propagator A. Then, the generating functionals W = InZ of the connected correlations are related by: 

W(A;j,M) = W(A;j) + i lnDct(A- 1 .A) = W(A;j) - ^Tr ln(A.A _1 ). (A.4) 

Next, using the fact that the integral of a derivative vanishes we obtain as in eas. l|51[) - (l52(l : 



, SW 1 
j — Sx + {M - A ).— — -S3 



SW SW S 2 W 



, SW SW S 2 W SW 
and -z-r-^r + t— = -— -. (A.5 
Sj Sj SjSj SM y ' 



5 3 Sj SjSj 

The first equation in (|A.5J| is simply the Schwinger-Dyson equation 1(52(1 while the second one is obtained by taking the 
derivatives of Z[j,M] defined in ea. l(A.l(l . Indeed, since Z[j,M] depends on two variables j,M, while the generating 
functional Z[j] defined in eg. 1(31(1 was a functional of the single test field j, we need two equations to fully define 



Z[j,M}. As pointed out in |Vasirev &: Kazanskii (1972) the second ea. l|A.5(l is actually a constraint equation which 
does not depend on the action S[<j>]. Next, we introduce the double Legendre transform T[tp, G] as in ea. 1)60(1 with 
N = 1. The functionals W and T are also related by the properties: 

Sip Sip Sj Sp sm , s 2 w ( s 2 r _ st _ s 2 r \ s 2 w _ s 2 r 

= — — • t. h - — ■ — — which yields 1 — „ „ . 



1 = T - = h ttt-"7 — which yields 1 = — — — ■ : 2— — 2 _^,_ .y + - . - - - .2 (A.6) 

and: 

„ Sip Sp Sj Sp> SM „ 5 2 W / <5 2 r „ <5 2 r \ (5 2 M^ „ <5 2 T 

= 77T = T7-^ + 7T7--F77 which yields = — -. — - - 2——.^ + — - .2— — (A.7) 



5G <5j <5G JM <5G SjSj \ScpSG SGSG'^J SjSM' SGSG 

where we used eq.(|HUl: M = 2ST/SG and j = ST/5p - 2{ST/5G).p. Combining both equations (|A~6)) - l|A~7)) we obtain: 

5 2 W ( S 2 T ST S 2 T f 5 2 T y 1 5 2 T \ 
~ SjSj'\S(pSip SG 8pSG\5G5GJ ' SGSp J ( 

Therefore, the equations of motion 1A.5|) also read: 

Of course, the first equation in ((A. 9(1 is consistent with the equation of motion we had derived in eq. J53J) for the first 
Legendre transform T[p], using the last property of 1(61(1 to relate T[p] to T[<p, G], and with ea. (|66|l written at the 
saddle-point. We can note that the first equation in (|A.9|I is easily solved as: T[<p, G] = S[<p] + (l/2)S s pG — F(G) 
where F(G) is an arbitrary functional of G. Then, F(G) is determined by the second equation l|A.9|) . Therefore, we 
can write T[tp, G] under the form: 

T[p,G] = S^ + ^TvlG^.G-l] -^Trm^^G) -f[G] with G^ 1 = = A" 1 + S 3 .p (A.10) 

This is indeed identical to the form 1(63(1 which we want to prove here for a cubic action. Thus, we have seen so far 
that the reduced effective action T only depends on G as T[ip, G] as written in ea. i(A.10(l is a solution of the first 
equation in l(A.9() . The reason why we explicitly subtract the terms Tr[A _1 .G — 1] and Trln(A _1 .G) from T is that 
for a Gaussian action (i.e. S3 = 0) we can easily solve for T[ip, G] which yields the expression i|A.10|) with T = 0, see 



also Vasil'ev & Kazanskii (1972) Then, substituting ea. l|A.10|) into the second equation (|A.9(1 we obtain: 
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Thus, this equation determines the reduced effective action f , up to a term f (°) which is independent of both <p and 
G. We can also check that T only depends on G and no longer on <p. This equation can be solved iteratively (see 
Vasil'ev & Kazanskii (1972) I using: 



G^G" 1 - 2 



S 2 f Y 



5GSG 



S 2 f S 2 f 8 2 f 

= GG + GG 2—— GG + GG 2^^ GG 2-rr^ GG 



5G5G 



5GSG 



SG6G 



(A.12) 



Graphically, this also reads: 

s 2 f y 1 

' SG5G ) 



G? _1 G _1 - 2- 



(A.13) 



where a line stands for G and a hatched circle for 8 2 T / 8G5G. This allows us to solve for T as a power-series over S3 
(except for the zero-term which is independent of G). Thus, we obtain for the first two terms: 



f« = _ 
12 



whence 



SG5G 



and f» = 1^3 



(A.14) 



Here the big dot with three legs is the vertex S3. We can note that these first few terms are 2PI vacuum diagrams 
with three-leg vertex £3 and propagator G. Moreover, they appear with the same multiplicity factors as in Wo{G), 
which is the term of order zero of the generating functional W(G;j) = X)^°=o V( n ')^™i™ corresponding to the action 



S[4>] 



.G- 



j|i-S'30 3 . In fact, these properties remain valid at all orders. First, it is easily seen from the recursion 



obtained from (|A. 1 lf> and (|A.13J| that all diagrams of f are 2PI vacuum diagrams, see |VasiTev fc Kazanskii (1972) 
(except for the possible zero-term T^°>). Indeed, the structure of the r.h.s. in eq. (|A.llfl is 2-line irreducible: 




(A.15) 



Second, we now need to show that f [G] contains all 2PI vacuum diagrams of Wq(G) with the same coefficients. From 
the definitions of T and T we have: 

f[G] = W\j, M] - j.<p - ±<p.M.<p - ^Tr[G.M] + S[<p] + ^Tr [G \G - l] - ^Trln (A" 1 ^) 



= W(A;j) - j.<p - \p.M. V - \tt[G.M] + S[<p] + ^Tr [G^ 1 .G - l] - ^Trln (a \g 
where we used ea. <|A.4|) in the second line. Besides, (J, M) are obtained from (cp, G) through: 



(A.16) 



M=2 o% = G » 1 - G - 1 - 2 m and i+^jj-* 



■A- 1 .<p + -S a [(fxp + G\. 



(A.17) 



Next, we can easily see from the definition of the Legendre transforms that Tg^Q = Tg 1= Q + S\.(p (this is obtained 
by noticing that the linear term S\.(j) of the action can be absorbed in the source term j.<p). Since this term S^.tp is 
included in the explicit term S[ip] in the r.h.s. in ea. (|A.10l) we see that T[G] is independent of both Si and ip. Therefore, 
we can evaluate the r.h.s. of ea. (|A.16p at the point <p = 0, Si = 0, which yields: 



f[G] = W(A;j) - \tt[G.M] + ^Tr .G - l] - iTrln (a" 
with: 



.G 



(A.18) 



] = -S 3 G and M 



A- 1 -G- 1 -2 — , hence A" 1 
oG 



G~ 



5T 



5T 



2— which yields A = G - 2G. — .G+ ...(A.19) 



SG 



SG 



The last equation in (|A.19|) is obtained by inverting the previous relation through an expansion over G. Substituting 



M and A into ea. (|A.18p we have 

6f 



T[G] = W(A;j) + Tr 
= W(K;j) + Tr 



G. 



SG 



iTrln l + 2G.g 



G^ 
G SG 



1 

2^ 



ST ST ST 

2G.— -2G. — .G.— 



(A.20) 
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where we developed the logarithm. First, we see from eas. (|A.19|) - (|A.20| that T[G = 0] = therefore the zero-term 
vanishes: — 0. Hence T only contains 2PI diagrams over G, as shown above from eas. (|A.ll() - i|A.15| ). so that 



r = r 



2PI 



Second, as noticed in Vasil'ev fc Kazanskii (1972)[ we can see that the first term of the series in ea. (|A.20|l 



cancels out with the second term in the r.h.s. of ea. (jA.20|l while higher-order terms are 2-particle reducible: 



Trln 1 + 2G 



ST 
SG 




(A.21) 



Here a line stands for G and a hatched circle for ST/SG. Therefore, we obtain: 



f [G] = f [G] 



2PI 



W[ A;j = ls 3 G 



2PI 



x ^ 

51 7T\ 



W„(A). -S 3 G 



(A.22) 



2PI 



where W^pj is the 2PI part of W (expressed in terms of G). In ea. ljA.22l) wc used the observation that if two functionals 
are equal their 2-particle irreducible parts (2PI) are also equal | |Vasil'ev fc K azanski i (1972)| ). Next, we note that for 
n > 1 the r.h.s. in eg. (|A.22[) contains at least one external leg attached to j which implies that the diagram is 1-particle 
reducible: 



n>l: W n .j n 



since j = an d A = G 



.(A.23) 



Here we singled out one branch attached to a weight j (shown by the cross). The big dot with three legs is again the 
vertex S3, the dashed line is the propagator A and the solid line is the propagator G. The large filled circle on the 
left is W„,.j™ _1 while the middle hatched circle in the right diagram is ST/SG. Therefore, only the first term Wq(A) 
in the series in ea. l|A.22|l (i.e. vacuum diagrams) contributes to the 2PI part of W . Next, we note that within these 
diagrams we may replace A by G since any insertion of a factor ST/SG yields a 2-particle reducible graph: 



2PI 



2PI 



2PI 



Therefore, we have shown (remember that W must be computed with Si = in ea. (|A.18|l and afterwards): 



f [G] = W (G) 



2PI 



2PI vacuum diagrams of the action S[ 



.G- 1 .<^ + -5; 



(A.24) 



(A.25) 



This completes the proof that for a cubic action S[4>] as given in ea. (|A.l() . the reduced effective action f introduced in 
eq. IjA.lOfl is given by the 2PI vacuum diagrams of the action S written in eq. (|A.25|) , with the same multiplicity factors 
as in Wo- 

For the case where we factorize an integer N in front of the action, as in (|31H . we simply need to make the 
following substitutions within the formulae derived in this appendix: S — » NS, W — > NW,j — > Nj, M — > NM, A -1 — > 
A^A -1 , T — » AT,(^ -> (p,G -> G/N and we obtain eqs. g^ - ^ . 
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